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Abstract 

The passage of Brownian motion through a bottleneck in a bounded domain is a 
rare event and as the bottleneck radius shrinks to zero the mean time for such passage 
increases indefinitely. Its calculation reveals the effect of geometry and smoothness 
on the flux through the bottleneck. We find new behavior of the narrow escape time 
. through bottlenecks in planar and spatial domains and on a surface. Some applications 

in cellular biology and neurobiology are discussed. 

cn- 1 Introduction 

' (-( '' The narrow escape problem is to calculate the mean first passage time (MFPT) of Brow- 

nian motion from a domain with mostly reflecting boundary to a small absorbing window. 
r-| ■ The MFPT, also known as the narrow escape time (NET), was calculated in [1]-[T2] for 

^ I small absorbing windows in a smooth reflecting boundary. Several more complex cases were 

considered in [6]- [8], such as the NET through a window at a corner or at a cusp in the 
boundary and the NET on Riemannian manifolds. The calculation of the NET in composite 
domains with long necks, as shown in Figure [H was attempted in [12] and [13] and ultimately 
. accomplished in [13] . The NET problem in a planar domain with an absorbing window at 

I the end of a funnel was considered in [15]. The case of planar domains that consist of large 



(N 



^ I compartments interconnected by funnel-shaped bottlenecks was also considered in [15 

■ In this paper we consider Brownian motion in two- and three-dimensional domains whose 

O , boundaries are smooth and reflecting, except for a small absorbing window at the end of a 

cusp-shaped funnel, as shown in Figure |2l The cusp can be formed by a partial block of a 
planar domain, as shown in Figure [31 The NET from this type of a domain was calculated 
in [15] only for the planar case. 

The results of [I]-[l2] for small absorbing windows in a smooth reflecting boundary of a 
domain Q can be summarized as follows. In the two-dimensional case considered in [7] the 
absorbing boundary dfla is a small window in the smooth boundary dfl that is otherwise 
reflecting to Brownian trajectories. The MFPT from x E Q to the absorbing boundary dQa, 
denoted fx-^ana^ is the NET from the domain Q to the small window dfla (of length a), 
such that e = 7r|(9r2a|/|9f2| = 7ia/\dQ\ <^ 1 (this corrects the definition in [7]). Because the 
singularity of Neumann's function in the plane is logarithmic the MFPT is given by 

T'x^dQa = "M; 1^ i*^^ +0(1) for a; G r2 outside a boundary layer near dVta. (1) 
nD 'n\o\La\ 



^Department of Applied Mathematics, Tel- Aviv University, Tel-Aviv 69978, Israel. Department of Math- 
ematics and Computational Biophysics, Ecole Normale Superieure, 46 rue d'Ulm 75005 Paris, France. This 
research is supported by an ERC-starting-Grant. 

^Department of Mathematics, Tel- Aviv University, Tel- Aviv 69978, Israel. 



1 





Figure 1 : A mathematical idealization of a cross section of smooth and sharp connections approx- 
imating the spine morphology: Left: The cross section is a composite domain that consists of a 
bulky head connected smoothly by an interface dO,i = AB to a narrow neck 0.2. The entire 
boundary is (reflecting), except for a small absorbing part d^a = CD. Right: A cross section 
of a sharp connection. 




Figure 2: Left: The planar (dimensional) domain Q,' is bounded by a large circular arc connected 
smoothly to a funnel formed by moving e apart two tangent circular arcs of radius Rc (i.e., AB = e). 
Right: Blowup of the cusp region. The red, green, and blue necks correspond to v± = 1, 0.4, and 
5 in (fT2|) . respectively. 
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Figure 3: Narrow straits formed by a partial block (solid disk) of the passage from the head to 
the neck of the domain enclosed by the black line. Inside the green circle the narrow straits can be 
approximated by the gap between adjacent circles. 



Surface of revolution formed when 

0.3+ 0.5 X x<2 

f(x) = \ ,0<x 

V 2.1 12500002 - {x-2.65f x < 4.103444186 

Olir ~ 4.103444186, is rotated about a horizontal axis. 

Figure 4: Narrow straits formed by a cone-shaped funnel. Axes are rotated 90° about the y-axis. 
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In the three-dimensional case the MFPT to a circular absorbing window d^la of small 
radius a centered at on the boundary dQ is given by 










AaD 


r L(o) + iv(o) 

271 


- a log a + o{a log a) 



(2) 



where 1/(0) and A^(0) are the principal curvatures of the boundary at the center of dVLa- 

However, the MFPT from a domain to an absorbing interface located at the end of a 
funnel, as shown in Figure O cannot be calculated by the methods of [6]- [8], [16], because the 
contribution of the singular part of Neumann's function to the MFPT in a composite domain 
with a funnel or another bottleneck is not necessarily dominant. Also the method of matched 
asymptotic expansions, used in [l]-[3], [17] for calculating the MFPT to the interface on a 
smooth boundary, requires major modifications for an interface at the end of a bottleneck, 
because the boundary layer problem does not reduce to the classical electrified disk problem 
[in]. Altogether different boundary or internal layers at absorbing windows located at the 
end of a cusp-like funnel are needed. The methods used in [12] and [12] for constructing the 
MFPT in composite domains of the type shown in Figure [Upright) are made precise here and 
the new method extends to domains of the type shown in Figure [T]( left). 

The new results of this paper are as follows. In Section |5] we prove that the MFPT to the 
narrow straits formed by a partial block of a planar domain (see Figures [2] and [3]) is given 

by 



r=W^^^^%^4^(l + «(!)) for.«|9fi|,i?„r„ (3) 

where Rc and Tc are the curvatures at the neck and e is the width of the straits. More general 
cases are also considered. In Section [3] we prove that the MFPT in the solid of revolution 
obtained by rotating the symmetric domain VL in Figure [2](left) about its axis of symmetry 
is given by 



3/2 



jgjd +»(!)) for a (4) 

where the radius of the cylindrical neck is a = e/2. In Section H] we consider Brownian 
motion on a surface of revolution generated by rotating the curve in Figure [2]( left) about its 
axis of symmetry. We use the representation of the generating curve 

y = r(x), A < X < 

where the x-axis is horizontal with x = A at the absorbing end AB. We assume that the 
parts of the curve that generate the funnel have the form 



r{x) = 0(y |x|) near x = 
(x - AV^" 

r(x) = a -\ z r^(l + o(l)) for z/ > near x = A, (5) 

^ ^ 1/(1 + ^ ^ ^ ^ 
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where a = \AB = e/2 is the radius of the gap, and the constant i has dimension of length. 
For V = 1 the parameter I is the radius of curvature Rc bX x = A. We prove that the MFPT 
from the head to the absorbing end AB is given by 



5(A) \{l + v)a 
~2D 



v/l+v 



^~^T7^ mf ' (6) 



sm 

1 + z/ 

where S is the entire unsealed area of the surface. In particular, for z/ = 1 we get the MFPT 

f (7) 

The case i/ = corresponds to an absorbing circular cap of small radius a on a closed surface. 
For a sphere we get the known result 

- 2i?% sinf 



sm^ 



where 9 is the angle between x and the south- north axis of the sphere and a = i? sin 5/2 (see 

mm)- 

If a right circular cylinder of a small radius a and length L is attached to the surface at 
z = A, the NET from the composite surface is given by 



5(A) V(l + y)a) " ^ SL ^ L 
sin 



2 



T = ' ^ , \ \ fora<£ (9) 



(see |T5j for a different derivation). We also find the NET and the exit probability when 
there are absorbing windows at the ends of narrow necks. These are related to the 
principal eigenvalue of the Laplacian in dumbbell-shaped domains that consists of heads 
interconnected by narrow necks, which, in turn, is related to the effective diffusion in such 
domains. In Section [5] we calculate the NET from composite domains that consist of a 
head connect by a funnel to a narrow cylindrical neck and calculate the principal eigenvalue 
composite and dumbbell-shaped domains. Finally, in Section [6] we calculate the mean time f 
for a Brownian needle to turn around in a tightly fitting planar strip. For a needle of length 
/o in a strip of width / it is given by 



Dr^/k{k-l)S Dr 



D 




(10) 



where Dr is the rotational diffusion coefficient and Dx is the translational diffusion coefficient 
along the needle. We close this article by providing several applications to cellular biology. 
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2 The MFPT to a bottleneck 



We consider the NET problem in an asymmetric planar domain, as in Figure E] or in an 
asymmetric version of the (dimensional) domain Q' in Figure [2J We use the (dimensional) 
representation of the boundary curves 

y' = r±{x'), A' < x' < for the upper and lower parts, respectively (11) 

where the x'-axis is horizontal with x' = A' at AB. We assume that the parts of the curve 
that generate the funnel have the form 



r±{x') = 0(a/|x'|) near x' = 
ix' - A'V'^"^ 

r±(x') = ±a' ± ^- + oil)) for > near x' = A', (12) 

^ ' z/±(l + i^±)i± ^ ^ =t ^ ^ 

where a' = \AB = e'/2 is the radius of the gap, and the constants £± have dimension of 
length. For v± = 1 the parameters l± are the radii of curvature at x' = A' . To simplify 
the conformal mapping, we first rotate the domain by 7r/2 clockwise to assume the shape in 
Figure [2](left). The rotated axes are renamed {x',y') as well. 

The NET of Brownian motion with diffusion coefficient D from a point x' = (x', y') inside 
the domain Q' with reflection at the boundary dQ' , except for an absorbing boundary dQ'^ 
at the bottom of the neck, is the solution of the boundary value problem 

DAu{x') = - 1 for X e (13) 

^^=0 ioT x' e on' - on', 

on 

u{x') =0 for x' e dQ'^. 

We convert to dimensionless variables by setting x' = i-^x, A' = £+A, the domain Q' is 
mapped into Q and we have (see (IT^ below) 

= el\Q\, \dQ'\ = i+\dQ\, \dQ'J = e' = i+\dQa\ = Ue- (14) 

Setting u{x') = u{x), we write (fT3l) as 

^Au{x) = - 1 for a; G (15) 



du{x) 



ioi xedVL- dVta 



dn 

uix) =0 for a; G dVta- 



2.1 Asymptotic analysis 

First, we consider the case v± = 1, = Re-, and /_ = r^, radius 1, and A has dimensionless 
radius rd Rc- This case can represent a partial block described in Figured With the scaling 
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Figure 5: The image in il^ = w{Q) of the necks AB in Figure [2] under the conformal mapping (jl6p . 
The funnel in Figure [2] (shown inverted in the third quadrant) is mapped onto the ring enclosed 
between the like-colored arcs and the large disk is mapped onto small red disk. The short black 
segment AB (of length e) is mapped onto the thick black segment AB (of length 2y/e + 0(e)). 



f lHj) the bounding circle B has dimensionless We construct an asymptotic solution for small 
gap e by first mapping the domain Q in Figure [Jt^left) conformally into its image under the 
Mobius transformation of the two bounding circles (thin line) into concentric circles. To this 
end we move the origin of the complex plane to the center of the right lower circle and set 

w = w[z) = — , (16) 

where 

2eRc + 2Rc + e^Rc + 2rc5 + 2rc 



a 



2{eRc + r, + R,) 

^ ^ye{8Rcrc + 4ei?2 + UeR^rc + 4:e^Rl + Sri + Ae'^R^rc + e^Rl + Aerl) 

2{eRc + rc + Rc) 



which maps the right lower circle into itself and VL is mapped onto the domain VL^ = w{Q) 
in Figure [5l The straits in Figure [2](left) are mapped onto the ring enclosed between the 
like-colored arcs and the large disk is mapped onto the small red disk. The radius of the 
small red disk and the elevation of its center above the real axis are 0{^/e). The short black 
segment of length e in Figure |2] is mapped onto a segment of length 2y/e + 0{e). 
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Setting u{z) = v{w) and e = 2rce/{Rc + r^), the system (fT5|) is converted to 

A,„v{w) = — — — r^TTT = = tor w e U,, 

D\w'{z)\^ 71) _ 1 + 0(5) |4 



dv(w) 

— - — =0 for G dQw — dQyj a 
on 

v{w) =0 for w e dQw,a- 

The MFPT is bounded above and below by that from the inverse image of a circular ring cut 
by lines through the origin, tangent to the red disk at polar angles 6 = Ci\/1 (brown) and 
9 = C2\/i (cyan) for some positive constants ci,C2, independent of e. Therefore the MFPT 
from r2 equals that from the inverse image of a ring cut by an intermediate angle 9 = c\/i 
(black) . 

The asymptotic analysis of f|T8l) begins with the observation that the solution of the 
boundary value problem f fTSj) is to leading order independent of the radial variable in polar 
coordinates w = re*^. Fixing r = 1, we impose the reflecting boundary condition at 6' = c\/l, 
where c = 0(1) is a constant independent of e to leading order, and the absorbing condition 
at ^ = TT. The outer solution, obtained by a regular expansion of f (e*^), is given by 

voie'') = Ai9-7r), (19) 

where A is yet an undetermined constant. It follows that 



-A. (20) 



d9 

To determine A, we integrate f lTSj) over the domain to obtain at the leading order 



89 
hence 



2ViA (21) 



A ~ (22) 



Now f|T9|) gives for 9 = c\/e the leading order approximation 



TT 



The following more explicit analysis was briefly summarized in [15] for the symmetric case 
v± = \,Rc = Tc and is explicitly given here for completeness. The leading order approxima- 
tion is obtained by an explicit integration of f llSp with respect to 9, 

'"^e f , f dri 



eve 
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so that 



V e 



D 



TX TX 



dip 



drj 



eve 

TT 



\eiv - 1 - e^'^Vil^ 



(vr — 77) dr] 



D J |e«') - 1 - e*'?yi|4' 



(25) 



First, we evaluate asymptotically the integral 

vr 



^2^ 



?7 (i?7 



?|4 



(26) 



by setting r/ = and noting that 



It follows that 



2sin2i# sin C VI 



+ 



0{CVi) for all r],e>0. (27) 



vr 



T] drj 



TT y e 



D J le^r, _ 1 _ e^r,y||4 D J \l + e + 0{ee)? D{c+l) 
- c 



CdC 



(1 + 0(v^)) . 



eve 



i2t 



Similarly, we obtain that 



TT 



drj 



dC 



c 



D J le^r, _ 1 _ eir,y||4 J \i + (^2^o{eC 

r: c 

eve 

where C = 0(1) is a constant, so that 



?/•2^|2 



= (l + 0(VI)), (29) 



To determine the value of the constant C, we note that implies that 



(30) 



dn 



dv 
86 



TT 



drj 



C 



D J - 1 - e*^yi|4 Dy/I V 
c\/l 



= fl + 0(v^)) (31) 
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and the integration of (1T8|) over gives 



2^1 



- dv {e'^) 



dn 



D 



(32) 



Now, dnH) and ([32]) imply that 4C = |fi|/2, so that the MFPT to the straits, f , is 

^^^11^^ + ^^^^^^^^^ + ^^^^^ fore«|9fiU+, (33) 
which is (!23|) . In dimensional units ( 133|) becomes 



2Zr^ '^'^ for £ < |,i?c,rc. 

In the symmetric case Rc = ( j34l) reduces to the result of [T5] 



(34) 



2DJe'/R, 



1 + o(l)) for < \dQ!lR,. (35) 



Next, we consider for simplicity the symmetric case i/+ = i/_ > 1, so i?c = = oo. After 
scaling the boundary value problem (fT3|l with (fT4|) . we can choose the bounding circles at 
A and B to have radius 1 and repeat the above analysis in the domain enclosed by the 
green curves, shown in Figure (jS]). The result ( 13^ becomes 



Tx\n'\ 



1 + o(l)) for e' < \dQ!\J+. (36) 



2Dv^ 

2.2 Exit from several bottlenecks 

In case of exit through any one of well-separated necks with dimensionless curvature 
parameters Ij and widths e^, we construct the outer solution ( |T9l) at any one of the N 
absorbing windows so that ( l20l) holds at each window. The integration of ( |T8l) over VL gives 
the following analog of (j2T|) . 



2.V^.^^ =-2l^V^.^~-^> (37) 
i=i j=i 



hence 



Equation (!33|) is then generalized to 



A ~ ^ . (3^ 



l + o(l)) for \dn\. (39) 
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Equations (1M|) - (136|) are generalized in a similar manner. 

To calculate the exit probability through any one of the necks, we apply the trans- 
formation (JT6ll separately for each bottleneck at the absorbing images dQw,ai , ■ ■ ■ , dflw,aN to 
obtain images Q^j for j = 1,2, . . . , N. Then the probability of exiting through dQw,a, is the 
solution of the mixed boundary value problem 



dv{w) 
dn 

v{w) 
v(w) 



for G flwi 
for u; G (9fi^, - |J dil^^ai 



(40) 



N 



1=1 



1 for G 9^]^„,a, 

for G dQ^^a,, j 7^ i- 



The outer solution, which is the exit probability through window dQ^^i, is an unknown 
constant pi. We construct boundary layers at each absorbing boundary dflw,aj for j i hj 
solving the boundary value problem in Qy^. , which is of the type shown in Figure [5] with a 
neck of width Sj. In each case the boundary layer is a linear function 



Vj {6) = 6ij — Aj {6 — re) for all j, 

such that 

t'j(O) ~ 5ij + AjTi = Pi for all j. 
To determine the value of the constant Pi, we note that 



dv (e^^) 



dn 
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-A 



(41) 



(42) 



(43) 



so the integration of (HUl) over gives 



N N 

v4j I dVtu,,aj I = ^ '^Aj ^JWj 



0. 



(44) 



The + 1 equations ( H2l) and ( l44l) for the unknowns pi, Ai, . . . , An give the exit probability 
from an interior point in the planar 



Pi 



. Je'/i. 



(45) 



3 The NET in a solid funnel-shaped domain 

We consider now the NET problem in the solid of revolution obtained by rotating the sym- 
metric domain Q' in Figure Et^left) about its axis of symmetry. The absorbing end of the 
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neck becomes a circular disk of radius a' = e'/2. Due to cylindrical symmetry of the mixed 
boundary value problem f llSp the MFPT in cylindrical coordinates centered on the axis of 
symmetry is independent of the angle. It follows that with the scaling f ll4p the boundary 
value problem f|T5|l in the scaled spatial domain Q can be written in cylindrical coordinates 
as 



Am 



1 du d'^u 
r dr dz'^ 



D 



(46) 



Equation f H6|) can be considered as a two-dimensional problem in the planar cross section 
by a plane through the axis of symmetry of Q in the (r, z) plane. Here r is the distance to 
the axis of symmetry of fl, the z axis is perpendicular to that axis and the origin is inside 
the cross section of Q, at the intersection of the axis with the tangent to the osculating circle 
to the cross section at the gap. Setting ui = ur^^"^, the MFPT equation f H^ takes the form 



dz"^ 



1/2 , 

D \ 4r2 



(47) 



in the cross section, with mixed Neumann- Dirichlet boundary conditions, as in the planar 
case. We assume that in dimensionless variables AB = e <^ 1 < so the funnel is a 

narrow passage. The transformation to the rotated and translated coordinates is given by 
f = r — 1 — e/2, z = —z + 1. Setting Mi(r, z) = u{r, z), equation f H7|) becomes 



d^u{f, z) 



dr 



+ 



d'^uir, z) 



dz"^ 



D 



1/2 



(4J 



3.1 Asymptotic solution 



We construct an asymptotic solution for small gap e by first mapping the cross section in 
the (r, 2;)-plane conformally into its image under the Mobius transformation flTB]) . 



(49) 



where a is given in (fT7|) for the symmetric case Rc = = 1. Setting u{C) = v{w), equation 
THl) becomes 



A„,v(w) 



D\w'{C)\^ 



( 



\ 



\ 



^ w \ a e 

7^e- + 1 + 77 

1 + aw 2 



1/2 



^ w ^ a e 
7^e- + 1 + 

l + aw 2 



(50) 



Because the normalized head of Figure |2]^left) is mapped into the narrow hot dog-shaped 
region in Figure |5] of width y/e at p = 1, we approximate 



w = e 



O(v^), 



w + a 
1 + aw 



l + 0(v^). 



(51) 
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We also have 



— 1 



k'(c)r 

so that (l47ll reduces to 



(1 + waf 



1 - a2 



^'2 



|1 - w + v^wl'* „ 

(1 + 0(7^)), (53) 



^,„._^i£(l^M)L(y^,i„|. (54) 

— w + Vewp\ 16 ' 



or equivalently, 



^" + ^—^rjv = ^- ^„ (1 + 0(v^)) . (55) 

Setting w = ^^(y — 16\/2)/D, we obtain the leading order equation 

y"iv) + ^ — ^^yiv) = 0. (56) 

The boundary conditions are 

y'{c^/I) = 0, 2/(7r) = 16V2. (57) 
The outer solution is the linear function 

2/outer(^) = Mr] + N, (58) 

where M and are yet undetermined constants. The absorbing boundary condition in ( !57l) 
gives 



2/outcr 



V) = Mvr + iV = 16^2. (59) 



A boundary layer correction is needed to satisfy the boundary conditions at the reflecting 
boundary at 77 = c^/e. To resolve the boundary layer at 77 = c^/e, we set 77 = ^/e^ and 
expand 



O(v^). 



Writing yhi{r]) = ^(Oj obtain to leading order the boundary layer equation 

+ ^^^l^2y no = 0, (60) 

which has two linearly independent solutions, 5^(0 ^^^d ^2(0 that are linear functions for 
sufficiently large C,- Initial conditions for Yi{C,) and ^2(0 can be chosen so that Y2{^) — )■ const 
as ^ ^ 00 (e.g., ^2(0) = -4.7, ^2(0) = -1, see Figure E]). Setting 
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Figure 6: Two linearly independent solutions of (|60p. The linearly growing solution Yi{^) (red) 
satisfies the initial conditions Yi{0) = 0, 1^(0) = 2. The asymptotically constant solution l2(0 
(blue) satisfies the initial conditions ^2(0) = —4.7, ¥2(0) = —1. The asymptotic value is ¥2(00) « 
-5. 



y^,(rj) = AY, [ + ByJ^ ] , (61) 



where A and B are constants to be determined, we seek a uniform approximation to y^t]) 
in the form yum{{v) = l/outcrl'?) + Vhiiv)- The matching condition is that AYi{ri/y/e) + 
+BYl{r|/^/e) remains bounded as ^ — 00, which implies A = 0. It follows that at the 
absorbing boundary 77 = vr we have 

VnnM =M7i + /3 -5B = I6V2 (62) 
At the reflecting boundary we have to leading order 

y'nnAcV^) = y'onUcV^) + VUcV^) = M + B^ = 0, (63) 



which gives 



The uniform approximation to v{w) is given by 



vunifipe'") = M ( - TT - , (65) 
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so that using (1521) . we obtain from (!65l) 



du 
dn 



drj 



w'iO 



2M 
= ^1 



1 + 0{V^)). 



(66) 



To determine the value of M, we integrate f|T5l) over Q, use fl66l) . and the fact that 



dS 



ire 



(67) 



to obtain M = —2i'^\Q\/D'iTe^^'^. Now (165!) gives the MFPT at any point x in the head as 



T = U[X) 



(68) 



The dimensional radius of the absorbing end of the funnel is a' = i+B/2 (see (fT4l) ). so ( l68|) 
can be written in physical units as 



r 



3/2 



V 

WD 



;i + o(l)) for a' < 



(69) 



where V = \VL'\ is the volume of the domain. 



3.2 Exit from several bottlenecks 



The generalization of ( l69l) to exit through N well-separated necks is found by noting that 
fl6711 becomes 



J. .7 = 1 



ana 

and the integration of ([13]) over f2' gives the compatibility condition (dimensional) 



Jan' dn' ^ 



(71) 



which determines 



^ = - nv^iv ' 3/2 - (72) 



Hence, using the dimensional a'j = ijej/2, we obtain 



^ = = 7^ ..3/2 - (73) 
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To calculate the exit probability from one of necks, we note that the boundary layer 
function is to leading order linear, as in Section 12.21 Therefore in the three-dimensional case 
the exit probability is given by 



Pi 



N 3/2 



EJV 



(74) 



4 Diffusion and NET on a surface of revolution 



We consider now Brownian motion on a surface of revolution generated by rotating the curve 
in Figure |2t^left) about its axis of symmetry and assume = = v and = = i. The 
projection of the Brownian motion from the surface to the z-axis gives rise to a drift. The 
backward Kolmogorov operator [22] of the projected motion, scaled with f lT^ . is given by 



C*u{z) = 1 , ,0, , 



u"{z) + 



r\z) 



r'{z)r"{z) 



r{z){l + r'{zY) {l + r'{zYy 
The operator C* corresponds to the Ito equation 

dz = a{z) dt + b{z) dw, 
where the drift a{z) and noise intensity b{z) are given by 



u'iz) 



(75) 



(76) 



r'{z) 



r'{zy\z) 



r{z){l+r'{zf) (1 + r'(^)^)2 j ' 



b{z) 



2D 



P{l+r"{z)) 



(77) 



and w{t) is standard Brownian motion on the line. The potential of the drift is A{z) = 
-Jla{t)dt. 

To calculate the MFPT from 2; = to the end of the funnel at 2; = A, we note that due to 
rotational symmetry the solution of the Andronov-Pontryagin-Vitt boundary value problem 
[22] for the MFPT u{z, 9) on the surface is independent of 9. Therefore the problem reduces 
to 



1 



d 



r{z)^Jl + r''\z)^z 



r{z) du{z) 



^l + r''\z) dz 

M (0) = m(A) = 0. 



£2 

D 



The MFPT is given by 



uiO) 



2ttD 



S{t) dt. 



(79) 
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Figure 7: The drift a{z) in ([77|) (left panel) and its potential A{z) (right panel) near the cusp. 
The projection of the Brownian motion on the axis of symmetry has an effective high barrier in 
the neck. 



where S{t) is the (scaled) area of the surface of revolution from z = t to z = 0, given by 



S{t) = 27ry" ris)^l + r'\s) ds. (80) 
t 

The main contribution to (I7U|) comes from A < t < A + 6 for a sufficiently small 5, such that 
6 ^ a (note that the singularity of l/r{z) near z = is integrable). Thus (1791) and (IHOl) give 
for u > 

u/l+u 

, = .(0) . TVTT^ . 5(A) ( (TT^ ' """" 



27rD J r{t) ^" 2D ■ 



1 + 



where iS = S{A) is the entire unsealed area of the surface. In particular, for = 1 we get 
the MFPT 

f (82) 

The case i/ = corresponds to an absorbing circular cap of a small radius a on a closed 
surface. For a sphere the solution of (1751) gives the known result 



2i?2 sinf 

Tx^an, = ^ log — ^, 83) 
u sm I 



where Q is the angle between x and the south-north axis of the sphere and a = i? sin 5/2 
(see [8]-[ll]). If a right circular cylinder of a small radius a and length L' = iL is attached 
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to the surface at 2 = A, then the integration in (IHT]) extends now to A — L, giving 



A A-L 



^5(A) ry/TTr^^ l;^ 

27rL) y r(t) 27iDa 2D' ^ ' 

A 

where the integral is given by fIST]) . flS^ . or f lS^ for the various values of u. Note that while 
f on the surface depends on the fractional power — 1^/(1 + 1^) of the neck's radius a, the power 
of a in the three-dimensional case is —3/2, as indicated in (169|) . 

The case = is not the limit of flHT]) . because the line (|5]) blows up. This case 
corresponds to a conical funnel with an absorbing circle of small radius a and length H (see 
Figure H]). We assume that the radius of the other base of the cone, b, is smaller than a, but 
that b <^ 5^/^. The generator of the cone is the line segment 

r{x) = a + C{x - L) for A - L < x < A, (85) 

where C is the (positive) slope. In this case f l51|l is replaced by 

„(o, ^ ^ / M . log ( 1 . 



2ttD J r{t) 2ttDC ^\ a 

(a + CL') log (^1 + + ^[(a + CL' 



A 

:i + c^) 



2DC^ 



'^2 - a"] 



which reduces to ( l84l) in the limit CL' <^ a and for a <^ CL' can be simplified to leading 
order to 

5(A) } Jl + r''\t) , 5(A)V1 + C2, CL' 

A 

(1 + C2)L'2 . CL' 



2D 



■log + 0(1). (86) 



Note that the last term in (186!) blows up as a — while that in ( 184|) does not. This is due 
to the degeneration of the NET problem in the cylinder, as noted in [6]. 



5 The principal eigenvalue in domains with bottlenecks 

The narrow escape time is related to the leading eigenvalues of the Neumann or mixed 
Neumann-Dirichlet problem for the Laplace equation in domains that consists of compart- 
ments and narrow necks. In domains that consists of compartments interconnected by narrow 
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Figure 8: A dumbbell-shaped domain consists of two large compartments 0,i and 1^3 connected by 
a narrow neck The bottleneck is the interval AB. 

necks the MFPT from one compartment to the other, as defined in [18], is to leading order 
(in the limit of shrinking neck) independent of the initial point of the escaping trajectory 
and is twice the MFPT from the compartment to the narrowest passage in the bottleneck 
(e.g., the interval AB in Figure [8]). Indeed, the reciprocal of this MFPT is to leading order 
the rate at which trajectories reach the bottleneck from the first compartment, so the recip- 
rocal of the MFPT is the lowest eigenvalue of the mixed Neumann-Dirichlet boundary value 
problem in the first compartment with Dirichlet conditions on the cross section of the neck. 

There is a spectral gap of order 1 from the smallest eigenvalue to the next one. It follows 
that long transition times of Brownian trajectories between compartments connected by 
bottlenecks are exponentially distributed and therefore the leading eigenvalues of Neumann's 
problem for the Laplace equation in a domain that consists of compartments interconnected 
by narrow necks are to leading order the eigenvalues of a Markov chain with transition rates 
that are the reciprocals of the MFPTs through the narrow necks, as is the case for diffusion 
in a potential landscape with several deep wells (high barriers) [201 EI] (see also [IS]). The 
evaluation of the leading eigenvalues of the Neumann problem for the Laplace equation in 
domains with bottlenecks reduces to the computation of the leading order eigenvalue for the 
mixed Neumann-Dirichlet boundary value problem for the Laplace equation in a domain 
with refiecting (Neumann) boundary except for a small absorbing (Dirichlet) window at the 
end of a funnel. Some estimates on the asymptotic behavior of the leading eigenvalue are 
given in [36], [37] and references therein. 

5.1 Eigenvalue of the mixed problem in domains with bottlenecks 

First we consider the principal eigenvalue of the mixed Neumann-Dirichlet problem for the 
Laplace equation in a composite domain that consists of a head Qi connected by a funnel 
to a narrow cylindrical neck The boundary of the domain is refiecting (Neumann) and 
only the end of the cylinder dQa is absorbing (Dirichlet). The left half of Figure |8] shows the 
composite domain and the absorbing boundary is the interval AB. In the three-dimensional 
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case the Dirichlet boundary d^la is a small absorbing disk at the end of the cylinder. The 
domain Qi is the one shown in Figure |2] and it is connected to the cylinder at an interface 
dQi, which in this case is the interval AB in Figure |21 It was shown in [T3] that the MFPT 
from X & Qi to dQa is given by 

fx^an. = fx^an. + ^ + (87) 

The principal eigenvalue of the mixed two- and three-dimensional Neumann-Dirichlet prob- 
lems in domains with small Dirichlet and large Neumann parts of a smooth boundary is 
asymptotically the reciprocal of the MFPT given in ( 1871) . Thus the principal eigenvalue Ai 
in a domain with a single bottleneck is given by 

1 

Ai ~ 



L2 l^^ilL 
Tx^an, + 777; + 



2D \dna\D 

where fx^ani is any one of the MFPTs given in ([I])-®, depending on the geometry of Qi. 

If a composite domain consists of a single head and TV well-separated bottlenecks of 
different radii and neck lengths, the derivation of ( 1731) shows that the reciprocal of the 
MFPT is the sum of the reciprocals of the NETs from a domain with a single bottleneck. 
That is, the principal eigenvalue Ap is given by 

N 

Ap~J]A,. (89) 

This can be interpreted as the fact that the total efflux is the sum of independent effluxes 
through the bottlenecks. 



5.2 The principal eigenvalue in dumbbell-shaped domains 

We consider now the principal eigenvalue of the Neumann problem in a dumbbell-shaped 
domain that consists of two compartments Qi and ^3 and a connecting neck ^2 that is 
effectively one-dimensional, such as shown in Figure El or in a similar domain with a long 
neck. We assume, as we may, that the stochastic separatrix (SS) in the neck is the cross 
section at its center. In the planar case it is the segment AB in Figure El This means that 
a Brownian trajectory that hits the SS is equally likely to reach either compartment before 
the other. Thus the mean time to traverse the neck from compartment fli to compartment 
^3 is asymptotically twice the MFPT fx-,ss from x e to the SS [IS]. This MFPT is to 
leading order independent of x E Qi and can be denoted tq^^^ss- 

First, we note that the mean residence time of a Brownian trajectory in Qi or in is 
much larger than that in Q2 when the neck is narrow. Second, we note that the first passage 
time Tx^ss for a; G r2i is exponentially distributed for long times and so is tx^ss for x G fl^ 
|22j . We can therefore coarse-grain the Brownian motion to a two-state Markov process (a 
telegraph process), which is in State I when the Brownian trajectory is in fli and is State 
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II when it is in ^l^. The state and the residence time there can be neglected relative to 
those in Qi and ^3. The transition rates from I to II and from II to I are, respectively, 



A 



i^ii 



2f, 



A 



2fn 



3^SS 



(90) 



These rates can be found from ( l88l) . with L half the length of the neck and SS = dQa- The 
radii of curvature i?c,i and -Rc.s at the two funnels may be different, and the domain is either 
Qi or ^3, as the case may be. The smallest positive eigenvalue A of the Neumann problem 
for the Laplace equation in the dumbbell is to leading order that of the two-state Markov 
process, which is A = —{Xj^u + Xn-^i) (see Appendix below). For example, if the solid 
dumbbell consists of two general heads connected smoothly to the neck by funnels (see ([1])), 
the two rates are given by 



1 



Xi^ii 
1 



A 



-V2 
-V2 



RcA 



3/2 



II^I 



R. 



c,3 



Rc,iD 



Rc,zD 



(1 + 0(1)) 



(1 + 0(1)) + 











4l) + 







(91) 
(92) 



Next, we consider the Neumann problem for the Laplace equation in a domain that consists 
of any number of heads interconnected by narrow necks. The Brownian motion can be 
coarse-grained into a Markovian random walk that jumps between the connected domains 
at exponentially distributed times with rates determined by the first passage times and exit 
probabilities, as described in Section 15.11 This random walk can in turn be approximated by 
an effective coarse-grained anisotropic diffusion, as done, for example, for atomic migration 
in crystals [3HI Ch.8, Sect. 2] and for effective diffusion on a surface with obstacles [T5] . 



6 A Brownian needle in dire straits 

As an application of the methodology described above, we study the planar diffusion of a 
stiff thin rod (needle) of length / in an infinite horizontal strip of width Iq > I. We assume 
that the rod is a long thin right circular cylinder with radius e <^ /q (Figure |9]). The planar 
motion of the rod is described by two coordinates of the centroid and the rotational angle 
9 between the axes of the strip and the rod. The y-coordinate of the center of the rod is 
measured from the axis of the strip. The motion of the rod is confined to the domain VL 
shown in Figure [Hb- The rod turns across the vertical position if it goes from the green to 
the blue domains or in the reverse direction. If 

^ = « 1, (93) 

the black window becomes narrow and the mean first passage times (MFPT) Tcreen-^Biack and 
Teiue^Biack, from the green or blue to the black segment, become much longer than those in 
the other directions. The former also become independent of the starting position outside a 
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boundary layer near the black segment. Thus the definition of the time to turn is independent 
of the choice of the green and blue domains, as long as they are well separated from the black 
segment. The neck near the black domain is the boundary layer region near 9 = 7r/2. We 
neglect henceforward the short times relative to the long ones. 

To turn across the vertical position the rod has to reach the black domain from the green 
one for the first time and then to reach the blue domain for the first time, having returned to 
the green domain any number of times prior to reaching the blue domain. Due to symmetry, 
a simple renewal argument shows that the mean time to turn, TBiue^creen, is asymptotically 
given by 

''"Blue^Green ~ 2TBlue^Black fo^^ j ^ 1- (9^) 

to 

The time to turn is invariant to translations along the strip (the x-axis), therefore it suffices 
to describe the rod movement by its angle 6 and the y coordinate of its center. The position 
of the rod is defined for 6 mod it. Therefore the motion of the rod in the invariant strip can 
be mapped into that in the {0,y) planar domain Q (see FigJOb): 

n=!^{9,y) : M<^^^^, 0<^<7r|. (95) 



6.1 The diffusion law of a Brownian needle in a planar strip 

In a rotating system of coordinates {X,Y,6), where the instantaneous X-axis is parallel to 
the long axis of the rod and the Y-axis is perpendicular to it, the diffusive motion of the rod 
is an anisotropic Brownian motion, and can be described by the stochastic equations 

X = \plD^wx 
Y = ^/2D^W2 
9 = v^«^3, 

where Dx is the longitudinal diffusion coefficient along the axis, Dy the transversal diffusion 
constant and Dr, the rotational diffusion coefficient. Due to the anisotropy, the rod makes in 
general larger excursions in the X-direction than in the F-direction and this usually charac- 
terized by the ratio — ^. In a fixed system of Cartesian coordinates (x, y), the translational 
Dx 

and rotational motion of the centroid {x{t),y{t)) and the angle of rotation 6{t) of the rod is 
governed by the Ito equations 

X = cos 

y = sm{9) \/2Dx Wi + cos(6') ^/2Dy W2 
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Figure 9: Top: Rod in strip . The strip width is and the rod length is / < ^o- The position of 
the rod is characterized by the angle 6 and the fixed coordinates x and y and the rotating system 
of coordinates {X,Y,9). Bottom: The motion of the rod is confined to the domain in the y,9 
plane. 



which can be put in the matrix form 



x{t) = B{9)w, 



where 





fx) 








y 















and 




B{9) = V2 



The probabihty density function of the rod in the product space ^2 x 
p{t, X, y, 9) dx = Pr{{x{t), y{t), 9{t)) ex + dx}, 



(96) 
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satisfies the Fokker-Planck equation 

dp{t, x) 



dt 



-V ■ J{t,x) 



where the flux is given by 



( [Dx cos' e + Dy sin' 0]^ + ^ [{Dx - Dy) sin 29] ^ ^ 



J{t,x) 



\ 



[Dx sin^ e + Dy cos' 0]^ + ^ [{Dx - Dy) sin 26] ^ 

dp 



Dr 



89 



(97) 



The boundary conditions are vr-periodic in 6, because the position of the rod is deflned 
modulo TT (note that J(t,x) is vr-periodic in 9). This means that the density p(t,x,y,9) is 
vr-periodic and the normal flux —Drdp{t,x,y,9)/d9 is vr-antiperiodic in 9. 



The MFPT Taiue^Biack is translation-invariant with respect to x and is, therefore, the 
solution u{9,y) of the boundary value problem 



D.. 



Dy{9) 



d'u{9,y) 



— 1 for 



(98) 



89' ' dy' 

where Dy[9) = Dx sin' 9 + Dy cos' 9 and f2i=f2n|0<— |, with the boundary conditions 



du 

dn 



"IT 



for [9, y) on the red boundary and at 6' = 
for \y\ < Iq — I, 



where the co-normal derivative of u{9,y) on the red boundary is given by 
du 



dn 



\/u{9, y) • n{9) for (6*, y) on the red boundary 



(99) 
(100) 

(101) 



and the co-normal vector n(^) is given by 

fi{9) = 



Dr 

Dy{9) 



n{9) 



(102) 



with n{9) - the unit normal vector to the red boundary. 
Introducing the dimensionless variables 

^, X Y xit) yit) 

Lq to h h 
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and the normalized diffusion coefficients 



n' - D' - n (0) - ^y^^^ 



we find tliat the domain ^2 in fl95|) is mapped into 

^^' = {(^,r/) : N<i^^^^^, 0<^<7r^. (103) 
To convert ( l98!) to canonical form, we introduce the variable 



m= I \l^d9', (104) 







which defines the inverse function 9 = 0{(p), and set u{6,y) = U^ip.rf) to obtain 

f/^^(y.,r/) + U,,{^,7^) = U^{^,v)^/D'r ^ - (105) 

The domain Q', defined in f llOSp . is mapped into the similar domain 

n" = /(^,„):|„l<i^(i^£)£i!lM, o<^<^wl (106) 



in the {ip, r]) plane. Because the the co-normal direction at the boundary becomes normal, 
so does the co-normal derivative. It follows that the no-fiux boundary condition fl99|) and 
the absorbing condition fllOOp become 

9U(jf, T]) _ ^ (QfipYr]) on d^l" (the red boundary in the scaled Figure (Hb) 
on 

= for 7] < - 

dip 2 

U[ip[^),v) = for |r/|<|, (107) 

respectively. The gap at ^ = 7r/2 is preserved and the (dimensionless) radius of curvature of 
the boundary at the gap is 



2n 



71 

v[2j _ 2Dx 



R' = = ^ — . (108) 



First, we simplify fll05p by setting 



9{^) = y^r^^^^l^, U{^,r^) = f{^)V{^,v) (109) 
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and choosing /((/?) such that f'{f) = \ f{ip)g{ip). Note that 

dD-^'\e) 



dd 



0. 



;iio) 



6l=0,7r/2,7r 



Equation fllOSp becomes 



1 



Next, we move the origin to the center of curvature of the lower boundary by setting 



111) 



C = -{v-R'-^)+^ 



if- if 



TT 



and use the conformal mapping (fT6|) . 



C-R'a 
R' -aC 



(112) 



with u = pe^^. We also have 

1 (1 + awY 



w'iO 



k'(C)P 



R' 



[l + waf 



\1 — w + \/ew 



AeR'^ 



-(1 + 0(7^)), 



(113) 
(114) 



The image of the domain VL is given in Figure [TO] and is similar to in Figure O except 





Rent' 



Figure 10: The image of the domam under the mapping (jll2|) . The values of the parameters 
are e = 0.01 with the approximation Dy ^ Dx- The domain is enclosed by the real segment AB 
and by the brown, black, and blue curves. The green and red curves are the images of arcs of the 
osculating circles at the narrow neck, as in Figure EJ 

for a small distortion near ijj = Cy/e, which we neglect, as we may. Setting V{ip, rf) = W{p, ifj), 
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fixing p = 1 in as in Section |2l and abbreviating W = W{ip, 1), equation f illip becomes 
to leading order 



where 



k{^) = f{^)D,{9{^))\,=,. (116) 



Using f|T8|) and neglecting terms of order 0{e), we rewrite (11151) as 

+ le^ni - V^) - 1|^ l^Hl - v^) - 1|^M^)' ^ ^ 

In view of (IllOp . tlie boundary conditions fll07p become 

WV(cv^) = 0, W{7r) = 0. (118) 
6.2 The asymptotic solution 

Tlie construction of the asymptotic expansion of the solution of the boundary layer equa- 
tion f lll7p is similar to that in Section 13.11 The outer solution of f lll7p is a linear function 
W^outer(V') = aip + b, where a and b are yet undetermined constants. The uniform approxima- 
tion is constructed as M/^uniform('^) = Wouteri'ip) + Wui'ip), where the boundary layer Wbi{ip) 
is a function Y{^) of the boundary layer variable C = "^Z V^- The boundary layer equation is 

which is simplified by the substitution Y{^) = Y{^) + l//i(0)/c(0) to 



^"(0 + ^!^^(0 = 0. (120) 



The boundary conditions (I118P become Y'{c) = and Y{oo) = l/h{0)k{0). The boundary 
layer equation f ll20p has two linearly independent solutions, Yi{^) and 1^(0 5 which are 
linear for sufficiently large ^. Initial conditions for Yi{^) and 5^2(0 can be chosen so that 
^2(0 ~^ const as ^ — > 00 (e.g., ^2(0) = —4.7, ^^'(O) = —1, see Figure[6]). Thus the boundary 
layer function is given by 

WuW = An (^^^ + BY2 (^^^ + C, (121) 

where A and B are constants to be determined and C is related to the constant l/h{0)k{0) 
and is also determined below from the boundary and matching conditions. 
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The matching condition is that Wbi(V') = ijP / \/^)+BY2 {jp / y/e)+C remains bounded 
as ^ — oo, which imphes A = Q. It follows that at the absorbing boundary ip = n we have 



W^unif (tt) =a7r + b' = 



(122) 



where the constant b' incorporates all remaining constants. At the reflecting boundary we 
have to leading order 



0, 



(123) 



which gives 



B 



ay/e 



b' 



The uniform approximation to W{u) is given by 



-air. 



Y'ic] 



(124) 



(125) 



so that using (1T09|) . f lUOj) . and (ITTSD . we obtain from (fT25|) 

'7r\\ dW{pe'^) 



du 
dn 



dip 



tl) = TT 



dip 
C=-i '89 



0=7r/2 



aM^,{l + 0{^e)). (126) 



Because W{uj) scales with 1/ f{ip) relative to V{ip, r]), we may choose at the outset f{ip{n/2)) 
1. 



Finally, to determine the value of a, we integrate fj98|) over fi, use (11261) . and the fact 



that 



dy = IqS, 



to obtain a = —{^lyR'/loDrV^- Now (I125P gives the MFPT at any point x in the head as 

7l\n\^/R 



E[t I x] = u{x) ~ W (pe'"^^ 



-an 



Reverting to the original dimensional variables, we get 

'TT 

n( 

E[r I x] 



;i + 0(v^) for e< 1. 



(127) 



TT I 1 

,2 



D 



X 



1 + 



/n-/ 



:i28) 



which is ([To 
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7 Discussion and conclusion 



This paper develops a boundary layer theory for the solution of the mixed Neumann-Dirichlet 
problem for the Poisson equation in geometries in which the methodologies of [1]-[I2] fail. 
These methodologies were used for the narrow escape problem. In the geometries considered 
here the small Dirichlet part is located at the end of narrow straits connected smoothly to 
the Neumann boundary of the domain. Additional problems related to Brownian motion in 
composite domains that contain a cylindrical narrow neck connected smoothly or sharply 
to the head are considered in [13]. These include the asymptotic evaluation of the NET, of 
the leading eigenvalue in dumbbell-shaped domains and domains with many heads intercon- 
nected by narrow necks, the escape probability through any one of several narrow necks, and 
more. 

Our results have applications in several areas. The first application is in neuroscience and 
concerns dendritic spines, which are believed to be the locus of postsynaptic transmission. 
Recognized more than 100 years ago by Ramon y Cajal, dendritic spines are small terminal 
protrusions on neuronal dendrites, and are the postsynaptic parts of excitatory synaptic 
connections. The spine consists of a relatively narrow cylindrical neck connected to a bulky 
head. The geometrical shape of a spine correlates with its physiological function [29]-|33]. 
Several physiological phenomena are regulated by diffusion in dendritic spines. For example, 
synaptic plasticity is induced by the transient increase of calcium concentration in the spine, 
which is regulated by spine geometry, by endogenous buffers, and by the number and rates 
of exchangers [311 |3ll |35l [231 [23] • Another significant function of the spine is the regulation 
of the number and type of receptors that contribute to the shaping of the synaptic current 
[2S]-[2S]- Indeed, the neurotransmitter receptors, such as AMPA and NMDA, whose motion 
on the spine surface is diffusion, mediate the glutamatergic-induced synaptic current. Thus 
dendritic spines regulate both two-dimensional motion of neurotransmitter receptors on its 
surface, and the three-dimensional diffusive motion of ions (e.g., calcium), molecules, proteins 
(e.g., mRNA), or small vesicles in the bulk. Our results give a quantitative measure of the 
effect of geometry on regulation of flux. 

In a biochemical context, the NET (eq. IBUjl accounts for the local geometry near an 
active binding site occluded by the molecular structure of the protein. This is the case of 
proton binding sites located on spike proteins, located on the viral envelope of the influenza 
virus and involved in membrane fusion [39]. Another application is that of the turnaround 
time of Brownian needle. Our result in Section [HI provides the precise time scale of the 
unraveling of a double strand DNA break confined between two-dimensional membranes 
[40j . The common feature of the geometries studied in this paper is the cusp-shaped narrow 
passage leading to the absorbing boundary. The main biological conclusion of our results is 
that this geometry is the main controller of the flux through biological narrow passages, an 
effect that is ubiquitous in biological systems. More specific applications of the NET from 
composite domains to dendritic spines are given in 
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8 Appendix 



The asymmetric random telegraph process jumps between two states, a and b, at inde- 
pendent exponentially distributed waiting times with rates Xa^b and Xb^a, respectively. 
The transition probability distribution function satisfies the linear differential equations (see 



http://en.wikipedia.org/wiki/Telegraph_process, [22] ) 



dP{a,t\x,to} I 1 , T-.r7 I T 

0^ = - Xa^bP{a, 1 1 X, to} + Xb^aP{b, 1 1 X, to} 

dP{b 1 1 X to} 

^^-^ = Xa^bP{a, 1 1 X, to} - Xb^aP{b, 1 1 X, to}, 

which can be written in the obvious matrix notation as p = Ap with 



(129) 



A 



'Xa^b Xb^a 
Xa^b —Xb^a 



The eigenvalues of A are with the normalized eigenvector (|, |)-^, and —{Xa^b + Xb^a) 
with the eigenvector (1,-1)-^. It follows that the nonzero eigenvalue of the system fll29p is 

X = Xa^b + Xb^a- 

Acknowledgment: The authors wish to thank F. Marchesoni for pointing out the factor 
1/2 in eq.§i). 
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